We need Python's plotting library, matplotlib. Your environment might load matplotlib automatically, but for this tutorial I'll load it explicitly using this convention. If you are unfamiliar with matplotlib, do the same as I do here, and everything that follows will work without modification.
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
We also might want to use scientific Python libraries. Finally, we'll import mr itself.
In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series # for convenience
import mr
Optionally, tweak the plotting styles.
In [29]:
plt.rcParams['figure.figsize'] = 10, 6
plt.rcParams['image.cmap'] = 'gray'
In [3]:
frames = mr.Video('../mr/tests/water/bulk-water.mov')
In [4]:
frames
Out[4]:
The Video constructor assumes brightfield video, inverting black and white. (See image below.) To process video of fluorescent particles, include the keyword argument invert=False.
In [5]:
mr.play(frames) # spawn a window playing the video
In [30]:
plt.imshow(frames[5], cmap=plt.cm.gray)
Out[30]:
Locate bright features with a diameter of about 11 px. The size must be an odd integer, and it is better to err on the large side. We will see why below.
In [32]:
f = mr.locate(frames[5], 11)
mr.locate returns a spreadsheet-like object called a DataFrame. It lists
(DataFrames are provided by the pandas module, and you can read more about them in the pandas documentation. They can easily be exported to a variety of formats.)
In [8]:
f.head() # shows the first few rows of data
Out[8]:
In [33]:
mr.annotate(f, frames[5])
Out[33]:
Many of these are spurious; they are fleeting peaks in brightness that aren't actually particles. There are many ways to distinguish real particles from spurrious ones. The most important way is to look at total brightness ("mass").
In [34]:
ax = f['mass'].hist()
ax.set_xlabel('mass')
ax.set_ylabel('count')
Out[34]:
In [11]:
f = mr.locate(frames[5], 11, minmass=2000)
In [35]:
mr.annotate(f, frames[5])
Out[35]:
There are more options in mr.locate, which you can review in the docstring using help(mr.locate) or mr.locate?.
In [13]:
help(mr.locate)
All functions in mr are documented in this way.
In [36]:
mr.subpx_bias(f)
Out[36]:
If we use a mask size that is too small, the histogram often shows a dip in the middle.
In [37]:
mr.subpx_bias(mr.locate(frames[5], 7))
Out[37]:
You can select:
frames[:]frames[100:200]frames[100:]frames[[100, 107, 113]]We'll locate features in the first 300 frames from this video. We use mr.batch, which calls mr.locate on each frame and collects the results.
In [16]:
f = mr.batch(frames[:300], 11, minmass=2000)
As each frame is analyzed, mr.batch reports the frame number and how many features were found. If this number runs unexpectedly low or high, you may wish to interrupt it and try different parameters.
If you have OpenCV installed (optional) you can view the circled features as a video.
In [17]:
mr.circle(f, frames[:300]) # opens a separate window
We have the locations of the particles in each frame. Next we'll associate each particle with an ID number and track it from frame to frame.
First, we must must specify a maximum displacement, the farthest a particle is likely to travel between frames. We must not underestimate how far a particle might travel, but we should choose the smallest reasonable value because a large value slows computation time considerably. In this case, 5 pixels is reasonable.
Second, we allow for the possibility that a particle might be missed in one frame and then seen again in the next frame. (Perhaps its "mass" slipped below our cutoff due to noise in the video.) Memory keeps track of disappeared particles and maintains their ID for up to some number of frames after their last appearance. We'll choose 3.
In [18]:
t = mr.link(f, 5, memory=3)
The result contains the information from the features f with an additional column, probe, identifying each feature with a label.
In [20]:
t.head()
Out[20]:
We have more filtering to do. Empheremeral trajectories -- seen only for a few frames -- are usually spurrious and never useful. The convenience function filter_stubs keeps only trajectories that last for a given number of frames.
In [23]:
t1 = mr.filter_stubs(t, 50)
# Compare the number of probes in the unfiltered and filtered data.
t['probe'].nunique(), t1['probe'].nunique()
Out[23]:
We can also filter trajectories by their appearance. At this stage, with trajectories linked, we can look at a feature's "average appearance" throughout its trajectory, giving a more accurate picture.
In [38]:
mr.mass_size(t1.groupby('probe').mean()) # convenience function -- just plots size vs. mass
Out[38]:
The probes with especially low mass or especially large size are probably out of focus or aggregated, respectively. It is best to experiment by trial and error, filtering out regions of mass-size space and looking at the results using mr.annotate and mr.circle. In the end, we need to separate the good probes from the spurrious ones, and it doesn't matter how we get it done.
In [47]:
condition = lambda x: (x['mass'].mean() > 2800) & (x['size'].mean() < 3.0) & (x['ecc'].mean() < 0.1)
t2 = mr.filter(t1, condition) # a wrapper for pandas' filter that works around a bug in v 0.12
In [48]:
mr.annotate(t2[t2['frame'] == 0], frames[0])
Out[48]:
Trace the trajectories.
In [49]:
mr.plot_traj(t1)
Out[49]:
In [50]:
d = mr.compute_drift(t1)
In [51]:
d.plot()
Out[51]:
In [52]:
tm = mr.subtract_drift(t1, d)
With the overall drifting motion subtracted out, we plot the trajectories again. We observe nice random walks.
In [53]:
mr.plot_traj(tm)
Out[53]:
Compute the mean squared displacement of each particle and plot MSD vs. lag time.
In [54]:
im = mr.imsd(tm, 100/285., 24) # microns per pixel = 100/285., frames per second = 24
In [62]:
im.plot(loglog=True, style='k-', alpha=0.1, legend=False) # black lines, semitransparent, no legend
plt.gca().set_ylabel(r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]');
In [57]:
em = mr.emsd(tm, 100/285., 24)
In [58]:
em.plot(loglog=True, style='o')
Out[58]:
We can easily fit this to a power law using a convenience function, fit_powerlaw, that performs a linear regression in log space.
In [63]:
em.plot(loglog=True, style='ro')
mr.utils.fit_powerlaw(em)
Out[63]:
In water, a viscous material, the expected power-law exponent $n = 1$. The value of $A = 4D$, where $D$ is the particles' diffusivity. $D$ is related to viscosity $\eta$, particle radius $a$, and temperature $T$ as
$D = \displaystyle\frac{kT}{6\pi\eta a}$.
For paritcles with a 1 $\mu\text{m}$ diameter in room-temperature water, $A \approx 1.74$. Our value above is not far off.